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With the Lattice Boltzmann method (using the BGK approximation) we investigate 
the dynamics of Hele-Shaw flow under conditions corresponding to various experimental 
systems. We discuss the onset of the instability (dispersion relation), the static proper- 
ties (characterization of the interface) and the dynamic properties (growth of the mixing 
zone) of simulated Hele-Shaw systems. We examine the role of reactive processes (be- 
tween the two fluids) and we show that they have a sharpening effect on the interface 
similar to the effect of surface tension. 

Viscous fingering occurs in the interfacial region between two fluids with dif- 
ferent viscosities when a highly viscous fluid confined between two plates with a 
narrow gap (Hele-Shaw geometry) is displaced by a fluid with relatively low vis- 
cosity. Most theoretical and numerical studies on Hele-Shaw systems start with 
Darcy's law Era assuming its validity under the simulation conditions, and thereby 
precluding discrimination between Navier-Stokes behavior or Darcy behavior of the 
system. With the Lattice Boltzmann (LB) approach U the discrimination problem 
can be addressed because the LB method starts with a kinetic theoretical analysis 
where the macroscopic description is not pre-established. 

While the LB method has now become a standard approach to investigate fln- 
gering processes, it has mostly been applied to the formation and growth of a single 
finger in a channel □. Fingering in large systems have been studied and analyzed 
with Lattice Gas Automata constructed to simulate porous media a u. Here we con- 
sider fingering in spatially extended Hele-Shaw systems. We investigate the effects 
of surface tension and reactivity between the two fluids, and we demonstrate their 
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role as determinant factors in the dynamics of the moving interface. The present 
paper gives a presentation of the essentials of the method and describes the main 
results of the work which will be discussed in detail elsewhere. 

We use the LB equatipn with the Bhatnagar-Gross-Krook collision operator 
(LBGK method; see e.g. a) in two-dimensional geometry emulating 3-D flows by 
introducing a Hele-Shaw drag term thereby simulating a system with a virtual cell 
gap in the third dimension. Using standard notation, the starting LBGK equation 
for multicomponent systems reads 

/r(r + c„t + l) - /f(r,t) = --^(/f - /;°^), (1) 

where a denotes the label of the component with density p'^{r) = fn" and 
velocity u'^(r) = -^jp^ ^ifi^i- Taking the long wavelength limit of the LBGK 
equation, one recovers the phenomenological equations of Navier-Stokes hydrody- 
namics (see e.g. 

Additional features can be incorporated in the LBGK equation to account for 
the physics of the problem to be considered. Here ^wc introduce (i) an interaction 
potential (initially used to simulate non-ideal fluids 0) for tunable miscibility ( "sur- 
face tension"); (ii) an external force for the Hele-Shaw drag term (a term that can 
also be used to incorporate tunable gravitational effect in density driven fingering), 
and (iii) a reactive term which quantifies changes in the "chemical" nature of the 
species which are susceptible to react with each other. 

(i) The interaction potential has the form V{r,r') = Gaa' {r,r') p'^ {r)p'^ (r') with 
C,{r,r') = C.{r-r') = G,,, for(r- r' ) = |r-r'| = c, and G,,,,- (r, r') = 
for |r — r'l > c, c being the modulus of the relative velocity between particles. The 
potential enters the dynamics through the modified local velocity 

p'^(r)u(r) = p'"(r)u(r) + T'^F'^(r,t), (2) 

with 

¥%r,t) = ~p%r) Y,G..' ^ p^'{r + c,)c.. (3) 

cr' i 

The amplitude (G) of the interaction term can be tuned from zero (miscible fluids) 
to practically total immiscibility. 

(ii) The external force (here the Hele-Shaw drag term) is introduced by a force term 
which modifies the momentum of species cr; the new momentum reads 

p%r) u-(r) = p^r) n^r) (1 - r") . (4) 

For the classical situation where fluid 1 invades fluid 2, the tunable drag coefficient 
is the ratio Pi/ (i2- 

(iii) Reactive processes are taken into account by a formulation similar to that used 
in reaction-diffusion equations (see e.g. □). Here we shall use a typical reactive term 
^ K pi P2 (pi ~ P2) for species 1 (and similarly for species 2 with 1 <^==> 2) where /« is 
the kinetic constant. This form restricts the reactivity effects to the mixing zoneB. 
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The LBGK model sketched above was used for binary fluids {a — 1,2) on 
the 2 — D square lattice with nine velocities oriented cdong the lattice axes and 
the bisecting directions plus a zero speed (D2Q9 model llJ). The simulations were 
performed on a L^; x Lj, = 1024 x 512 lattice with the following initial conditions: 
fluid 1 with lowest viscosity is injected from the upper boundary of the system 
filled with the highly viscous fluid 2; periodic boundary conditions are imposed on 
the vertical boundaries; the initial planar interface is perturbed with white noise 
to trigger the instability. As time evolves, the highly viscous fluid 2 is displaced by 
fluid 1. The color code used in Fig.l is from dark gray to light gray for p2 = — > 1 
so that the pictures show "dark" fluid invading from above the system initially flUed 
with "hght" fluid. 

Preliminary simulations were performed with a single phase fluid subject to a 
density gradient (Vp) imposed along the vertical axis (y), and the average velocity 
Vy was measured for various values of the gradient: Vy was found to be a linear 
function of the pressure gradient (oc Vp) according to Darcy's law. 

The top panel in Fig.l illustrates the case of miscible fluids (G = in Eq.(|3[)) 
where a distinguishable mixing layer can be quantified to measure the mixing 
length Lmix as a function of time t for increasing values of the Peclet number 
(Pe = LxVy/D where D is the diffusion coefficient). In agreement with experimen- 
tal observation and numerical simulation EJ, we observe a transition from the short 
time diffusive regime where Lmix scales as with Q!~0.5toa~l behavior 
characteristic of the non-linear regime. 

The middle panel of Fig.l shows the effect of surface tension on the finger topol- 
ogy (compare with upper panel). Increasing the amplitude (G) of the interaction 
term reduces drastically the width of the mixing layer which eventually vanishes, 
and the finger wavelength is modified correspondingly. A quantification of this effect 
is given in Fig. 2 where we show typical dispersion curves. 

Reactive processes between two moving miscible fluids are expected to modify 
the fingering mechanisms as shown by phenomenological approaches based on a 
reaction-diffusion scheme □, and experimental studies have indeed demonstrated 
that the interface properties depend strongly on chemical reactions t3. We have 
implemented a reactive mechanism by incorporating in the LB scheme a reaction 
term as described in (iii) above. The resulting effects are illustrated in the lower 
panel of Fig.l which shows (i) a sharpening of the intcrfacial zone, and (ii) a mod- 
ification of the topology of the fingering (compare with top panel) as quantified by 
the dispersion curves shown in Fig. 2. 

The purpose of this short paper is to demonstrate the applicability of the lattice 
Boltzmann method to the investigation of viscous fingering phenomena in spatially 
extended systems. In particular we have shown how surface tension and reactive 
processes can modify the dynamics and the characteristics of the fingering pattern. 
A quantitative analysis will be presented in a forthcoming paper with a detailed 
description of the work. 
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Fig. 1. LB simulations of viscous fingering (after t = 6000 time steps) in miscible fluids, G = 
(top), fluids with surface tension, G = 0.15 (middle), and miscible (G = 0) reactive fluids (bottom). 
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Fig. 2. Dispersion curves for (a) miscible fluids (G = 0; dots), (b) fluids with surface tension 
(G = .05; circles), and (c) miscible reactive fluids (k = .001; squares) showing the effect of surface 
tension and of reactive processes on the most unstable wavelength and on the cutoff wave number. 
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